How large is 1ha in size? Below is an example of a 1ha plot at Moore Reef. Note that this is only for demonstration purposes:

# data packages
library(tidyverse)
library(units)
library(foreach)

# spatial packages
library(sf)

# mapping packages
library(tmap)
library(leaflet)

set_restoration_plot_centres <- function(input = NULL, width = NULL, length = NULL) {

  # Calculate the coordinates of the rectangular polygon
  x <- sf::st_coordinates(input)[1, 1]
  y <- sf::st_coordinates(input)[1, 2]

  # set parameters
  x_min <- x - (width / 2)
  x_max <- x + (width / 2)
  y_min <- y - (length / 2)
  y_max <- y + (length / 2)

  polygon <- sf::st_polygon(list(rbind(c(x_min, y_min), c(x_min, y_max), c(x_max, y_max), c(x_max, y_min), c(x_min, y_min)))) |>
    sf::st_sfc(crs = 20353)

  return(polygon)
}

habitats <- c("Plateau", "Back Reef Slope", "Reef Slope", "Sheltered Reef Slope", "Inner Reef Flat", "Outer Reef Flat", "Reef Crest")

moore_geomorphic <- st_read("/Users/rof011/GBR-coral-restoration/data/Moore-Reef-20230906030257/Geomorphic-Map/geomorphic.geojson", quiet=TRUE) %>%
   sf::st_transform(crs = sf::st_crs("EPSG:20353")) %>%
   dplyr::group_by(class) %>%
   dplyr::mutate(habitat_id = paste0(
   gsub(" ", "_", class), "_",
   sprintf(paste0("%0", ceiling(log10(max(1:length(class)))), "d"), 1:length(class)))) %>%
   sf::st_make_valid() %>% 
#  mutate(habitat_id=as.factor(habitat_id)) %>%
  dplyr::group_by(habitat_id, class) %>% 
  summarise(geometry = st_union(geometry)) %>% 
  mutate(area = round(st_area(geometry))) %>%
  filter(class %in% habitats) %>%
  drop_na(class) 
  

moore_plots_centroids <- moore_geomorphic %>%
  filter(area > set_units(10000,"m^2")) %>%
  filter(class %in% c("Reef Slope", "Sheltered Reef Slope", "Back Reef Slope")) %>% 
  drop_na(class) %>%
  st_centroid() %>% 
  st_cast("POINT")

moore_plot_outputs <- foreach(i=1:nrow(moore_plots_centroids), .combine="c") %do% {
  
  tmp <- set_restoration_plot_centres(moore_plots_centroids$geometry[i], width=100, length=100)
  
  tmp
  
}


generate_random_points <- function(polygon, n) {
  st_sample(polygon, size = n)
}


moore_plot_outputs3 <- moore_plot_outputs |> st_as_sf() |>  slice_sample(n=1)

moore_plots_centroids_small <- st_as_sf(st_centroid(moore_plot_outputs3))

#moore_plot_small <- set_restoration_plot_centres(moore_plot_outputs3, width=10, length=10) #|> st_as_sf() |> mutate(id=moore_plot_outputs3$id)

  width=10
  length=10
 # Calculate the coordinates of the rectangular polygon
  x <- sf::st_coordinates(moore_plots_centroids_small)[1] |> as.numeric()
  y <- sf::st_coordinates(moore_plots_centroids_small)[2] |> as.numeric()

  # set parameters
  x_min <- x - (width / 2)
  x_max <- x + (width / 2)
  y_min <- y - (length / 2)
  y_max <- y + (length / 2)

  moore_plot_small <- sf::st_polygon(list(rbind(c(x_min, y_min), c(x_min, y_max), c(x_max, y_max), c(x_max, y_min), c(x_min, y_min)))) |>
    sf::st_sfc(crs = 20353) 

all_random_points <- st_sample(moore_plot_small, size = 5000) |> st_as_sf(crs=20353)
   
coral_locations <- all_random_points %>% 
  mutate(coral_id=seq(1,n(),1)) %>%
  mutate(coral_id=as.factor(coral_id)) %>%
  st_buffer(dist = 0.01)

coral_locations_yr5 <- all_random_points %>% 
  slice_sample(n=1000) %>%
  mutate(coral_id=seq(1,n(),1)) %>%
  #mutate(radius=2.5) |> 
  st_buffer(dist = 0.01+(0.025*4)) 
 

#moore_satellite_inset <- raster::crop(moore_satellite, st_bbox(moore_plot_outputs2 |> st_transform(3857)))

# make an internal box of 10x10 to speed things up

map <- tm_shape(moore_plot_outputs3, is.master=TRUE) +
  tm_borders(col=NA) +
tm_tiles("Esri.WorldImagery", group = "<b> [Seascape]</b> satellite map", alpha = 0.6) +
tm_shape(moore_satellite_inset, raster.downsample=FALSE, alpha=0.6) +
  tm_rgb(r=1, g=2, b=3) +
tm_shape(moore_geomorphic, id="area", name = "<b> [Seascape]</b> habitats") +
  tm_borders(col = "black", lwd = 0.2) +
  tm_fill("class", name = "Benthic habitats", id="area", alpha = 0.6) +
 tm_shape(moore_plot_outputs3, name = "<b> [Restoration]</b> 1ha plot") +
   tm_borders(col="darkred", lwd=1.5) +
 tm_shape(moore_plot_small, name = "<b> [Restoration]</b> 100m2 plot", is.master=TRUE) +
   tm_borders(col="darkred", lwd=1.5) +

  tm_shape(coral_locations, id="coral_id", legend.show=FALSE,   name = "<b> [Corals]</b> Year 1") +
   tm_polygons(col="darkblue", alpha=0.4) +
 tm_shape(coral_locations_yr5, id="coral_id", legend.show=FALSE, name = "<b> [Corals]</b> Year 5") +
   tm_polygons(col="aquamarine3", alpha=0.4) +

 tm_view(set.zoom.limits=c(5,100), bbox=st_bbox(moore_plot_small))

  map |> tmap::tmap_leaflet() |>
    leaflet::addProviderTiles('Esri.WorldImagery', group = "<b> [Seascape]</b> satellite map", options=leaflet::providerTileOptions(maxNativeZoom=18,maxZoom=100)) |>
    leaflet::addProviderTiles('Esri.WorldTopoMap',  group = "<b> [Seascape]</b> base map", options=leaflet::providerTileOptions(maxNativeZoom=19,maxZoom=100)) |>
    leaflet::addLayersControl(position="bottomleft", overlayGroups=c(
      "<b> [Seascape]</b> satellite map", "<b> [Seascape]</b> habitats", "<b> [Restoration]</b> plots", 
      "<b> [Restoration]</b> 1ha plot", "<b> [Restoration]</b> inset plot", "<b> [Corals]</b> Year 1", "<b> [Corals]</b> Year 5"),
                              options=leaflet::layersControlOptions(collapsed = FALSE)) |>
    leaflet.extras::addFullscreenControl(position = "bottomleft", pseudoFullscreen = FALSE)